#Landscape Connectivity and burrowing owl artificial burrow system locations

Research introduction and areas of interest

My research focuses primarily on pedigree analyses to help understand population dynamics in a burrowing owl population that has been monitored for approximately 20 years. Thus, my research will combine field work that comprises the trapping, banding, and and obtaining different physiological metrics of a wild burrowing owl population during the next couple of breeding seasons with pedigree analysis of a long-term dataset to understand aspects of kin structuring and reproductive success within that population.

I am interested in how kin relationships are structured spatially and temporally in a populationm, quantifying lifetime reproductve success and examining if and how production of offspring and recruitement into the breeding population differs by landscape type, and what are the patterns in burrow occupancy as it relates to nestling production and recruitement.
## Problem to be solved The study area is located in the Morley Nelson Snake River Birds of Prey National Conservation Area and within the study area there have been more than 350 different artificial burrow systems (ABS) constructed and installed throughout the landscape with the intent that the burrowing owl population that migrate to the area during the breeding system will utilize the ABS as their breeding site. The majority of the study area is comprised of sage steppe and agriculture with rolling hills. As human peturbation (agriculture, recreation, etc.) and encroachment of invasive plant species increases potentially limiting future habitat, how is the connectivity between ABS affected.

Psuedocode

Steps lining out how the problem will be addressed

#Analysis/processing step 1: 
##Provide a map detailing the Burrowing Owl population study area in southwestern Idaho
###-create raster of Idaho
###-create polygon's of Idaho and Ada and Elmore counties
###-crop raster to only include Ada and Elmore counties
###-insert geographic points of all burrowing owl locations
###-combine all of the above into a single map

#Analysis/processing step 2:
##Calculate the least cost distances (LCD) between every ABS location with the package 'gdistance'
###-create an empty raster of Ada and Elmore counties and set all cells within that raster so there is unity between each cell (cell value of 1) Do this to eliminate bias
###-insert ABS coordinates 
###-calculate LCD of ABS using the 'costDistance' function
###-Analyze output

#Analysis/processing step 3:
##Calculate the least cost distances (LCD) between every ABS location with the package 'igraph'
###-create a hexagonal grid of Ada and Elmore counties
###-create a 'igraph' network of cells based on the hexagonal grid
###-insert ABS coordinates 
###-calculate LCD of ABS using the 'distances' function
###-Analyze output

#Profile and Benchmark the two different r packages used 'gdistnace/igraph'

Install and Load required packages

Introduce the study area (State of Idaho)

###
#Create the state of Idaho and elevation Map
###

#Get the Global Administration boundries with state lines of the United States
USA_0 <- getData(name = 'GADM', country = 'USA', level = 1)
#plot(USA_0)

#create the state of Idaho
idaho <- c("Idaho")

#extract the state of idaho from the GADM
ID = USA_0[match(toupper(idaho),toupper(USA_0$NAME_1)),]
#plot(ID)

#crop the GADM of US to only indclude the state of idaho
idaho_0 <- crop(USA_0, ID)
#plot(idaho_0)

#Get the SRTM with 90m resolution elevation data of the area around Idaho using longitude and latitude
elevID_1 <- getData("SRTM", lon = -115, lat = 44)
elevID_2 <- getData('SRTM', lon = -117, lat = 46)
elevID_3 <- getData("SRTM", lon = -118, lat = 43)
elevID_4 <- getData("SRTM", lon = -113, lat = 47)

#mosaic the SRTM's so that all of Idaho is covered with elevation raster
elevID_mosaic <- mosaic(elevID_1, elevID_2, elevID_3, elevID_4, fun = mean)
#plot(elevID_mosaic)
#plot(idaho_0, add = TRUE)

#crop the mosaic of SRTM data with the boundry polygon of Idaho and plot new map
elevIdaho <- crop(elevID_mosaic, idaho_0)

###add counties where 'ABS' burrows are located
#read in idaho counties shape file
idaho_counties <- readOGR("~/library/Mobile Documents/com~apple~CloudDocs/BSU/Graduate - MS Raptor biology/r_Projects/HES598/datatemp/original/idaho_cty.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/brentclark/Library/Mobile Documents/com~apple~CloudDocs/BSU/Graduate - MS Raptor biology/R_Projects/HES598/datatemp/original/idaho_cty.shp", layer: "idaho_cty"
## with 44 features
## It has 17 fields
#create separate map of Ada & Elmore counties
counties_AE <- subset(idaho_counties, NAME == "Ada" | NAME == "Elmore")

#plot map with elevations, Idaho, and counties of interest
plot(elevIdaho, main = "Idaho Elevation (SRTM 90m Res)", xlab = "Longitude", ylab = "Latitude", legend.args = list(text = 'Elevation (m)', side = 4, line = 3, font = 2))
plot(idaho_0, add = TRUE)
plot(counties_AE, add = TRUE)

Introduce study area at the county level (Ada & Elmore Counties)

###
#create an elvation map of Ada & Elmore counties 
###

#crop Idaho raster with county pologons to create only counties raster
elevCounties <- crop(elevIdaho, counties_AE)

#plot map of Ada & Elmore counties and elevations
plot(elevCounties, main = "Ada & Elmore County Elevation (SRTM 90m Res)", xlab = "Longitude", ylab = "Latitude", legend.args = list(text = 'Elevation (m)', side = 4, line = 3, font = 2))
plot(counties_AE, add = TRUE)

Create each of the artificial burrow locations based on their geographic coordinates

###
#create the artificial burrows as a spatial object
###

#read excel data on artificial burrow locations (Lat & Long) into new data_frame
artBurrow_latlong <- read_excel("~/library/Mobile Documents/com~apple~CloudDocs/BSU/Graduate - MS Raptor biology/r_Projects/HES598/datatemp/original/Belthoff Burrows - latlong_revised.xlsx")

#from artificial burrow location data_frame create spatial table 
#artBurrow_spatialTable <- artBurrow_utm %>% st_as_sf(coords = c("EASTING", "NORTHING"), crs = 4326)
artBurrow_spatialTable <- artBurrow_latlong %>% st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326)

#creates burrow location points based on 'artBurrow_spatialTable'
artBurrow <- st_geometry(artBurrow_spatialTable) 
artBurrow <- st_set_crs(artBurrow, 4326)
plot(artBurrow, main = "ABS Geographic Loacations", cex = 1.5, col = "black")

Overlay the artificial burrow locations onto study area

###
#Overlay 'ABS' locations onto Ada & Elmore county Elevation map ('elevCounties')
###

#plot map of Ada & Elmore counties and elevations with 'ABS' locals
plot(elevCounties, main = "'ABS' locations - Ada & Elmore County Elevation (SRTM 90m Res)", xlab = "Longitude", ylab = "Latitude", legend.args = list(text = 'Elevation (m)', side = 4, line = 3, font = 2))
plot(counties_AE, add = TRUE)
plot(artBurrow, cex = 0.8, col = "black", add = TRUE)

Create a spatial points data frame of the ABS geographic locations

#remove 'ABS_NAME' column from 'artBurrow_latlong' data frame
artBurrow_latlong_noNames <- subset(artBurrow_latlong, select = -c(ABS_NAME))

#create object for coordinates
longLat <- artBurrow_latlong_noNames[,c(1,2)]

#createa a SpatialPointsDataframe()
artBurrow_spdf <- SpatialPointsDataFrame(coords = longLat, data = artBurrow_latlong, proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 "))

#find duplicated points in SpatialPointsDataFrame
dupPoints <- zerodist(artBurrow_spdf, zero = 0.0)
#remove duplicated points in SpatialPointsDataFrame
artBurrow_spdf <- remove.duplicates(artBurrow_spdf, zero = 0.0)

Spatial continuity of distances between Artificial Burrow Systems ‘ABS’

Least Cost Distance (LCD) Calculation with two packages will be used and compared:

1. ‘gdistance’ package

2. ‘igraph’ package

‘gdistance’ package using the function ‘costDistance’

Create empty raster from counties shape file

#read in Idaho counties shape file
counties <- st_read("~/library/Mobile Documents/com~apple~CloudDocs/BSU/Graduate - MS Raptor biology/r_Projects/HES598/datatemp/original/idaho_cty.shp") %>% st_transform(crs = 32619)
## Reading layer `idaho_cty' from data source `/Users/brentclark/Library/Mobile Documents/com~apple~CloudDocs/BSU/Graduate - MS Raptor biology/R_Projects/HES598/datatemp/original/idaho_cty.shp' using driver `ESRI Shapefile'
## Simple feature collection with 44 features and 17 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -117.243 ymin: 41.98818 xmax: -111.0435 ymax: 49.00091
## epsg (SRID):    NA
## proj4string:    +proj=longlat +ellps=GRS80 +no_defs
#isolate Ada & Elmore countie pologons from Idaho county shape file
studyCounties <- subset(counties, NAME == "Ada" | NAME == "Elmore")

#bring in spatial coordiantes of 'ABS' locals
coords <- artBurrow_spdf

#set 'ABS' spatial coordinates from lat long to UTM 19N
coords <- st_as_sf(coords, coords = c("Longitude", "Latitude"), crs = 4326, stringAsFactors = FALSE) %>% st_transform(crs = 32619)

#create raster of 'studyCounties' spatial object
studyCounties_ras <- raster(x = extent(studyCounties), nrow = 250, ncol = 250)

#assign all cells to a value of 1 within the raster
studyCounties_ras <- rasterize(x = studyCounties, y = studyCounties_ras, field = 1)

Convert ‘studyCounties_ras’ raster into a TransitionLayer with ‘gdistance’ package

#create an object of the Transition class
studyCounties_ras_tr <- transition(studyCounties_ras, transitionFunction = mean, 8)

#Geo correct transition object
studyCounties_ras_tr <- geoCorrection(studyCounties_ras_tr, type = "c")

Least Cost Distance calculation with ‘gdistance’

#Conduct least cost distance calculation on all 'ABS' locals
studyCounties_leastCost <- costDistance(studyCounties_ras_tr, fromCoords = as(as_Spatial(coords[1,]), "SpatialPoints"), toCoords = as(as_Spatial(coords[-1,]), "SpatialPoints"))
#put LCD calculation into data.frame
studyCounties_leastCost_df <- data.frame(studyCounties_leastCost, row.names = "Dist. from Initial Point (m):")
head(studyCounties_leastCost_df)
##                               X1 X2 X3       X4       X5       X6       X7
## Dist. from Initial Point (m):  0  0  0 772.6167 772.6167 772.6167 772.6167
##                                     X8       X9      X10      X11      X12
## Dist. from Initial Point (m): 1024.549 4740.604 4740.604 4740.604 4740.604
##                                    X13      X14      X15      X16      X17
## Dist. from Initial Point (m): 4740.604 4740.604 5413.481 56064.98 55040.43
##                                    X18     X19      X20      X21      X22
## Dist. from Initial Point (m): 55040.43 54788.5 39319.34 39319.34 39319.34
##                                    X23      X24      X25      X26      X27
## Dist. from Initial Point (m): 40091.96 39319.34 40091.96 40091.96 39319.34
##                                    X28      X29      X30      X31      X32
## Dist. from Initial Point (m): 52487.47 52487.47 52487.47 52487.47 52487.47
##                                    X33      X34      X35      X36      X37
## Dist. from Initial Point (m): 51462.92 51462.92 51462.92 51462.92 19533.71
##                                    X38      X39      X40      X41      X42
## Dist. from Initial Point (m): 19533.71 52991.33 53243.27 52991.33 52470.65
##                                    X43      X44      X45      X46      X47
## Dist. from Initial Point (m): 53243.27 52991.33 52470.65 52470.65 52470.65
##                                    X48      X49      X50      X51      X52
## Dist. from Initial Point (m): 53243.27 49380.18 52991.33 54015.88 53763.95
##                                    X53      X54     X55     X56     X57
## Dist. from Initial Point (m): 51698.03 51698.03 50152.8 50152.8 50152.8
##                                    X58      X59     X60      X61      X62
## Dist. from Initial Point (m): 49380.18 49380.18 48859.5 48086.88 48086.88
##                                    X63      X64      X65      X66      X67
## Dist. from Initial Point (m): 47834.95 28250.78 28250.78 28250.78 28502.71
##                                    X68      X69      X70      X71      X72
## Dist. from Initial Point (m): 19533.71 19785.64 13638.35 60969.43 60196.82
##                                    X73      X74      X75     X76     X77
## Dist. from Initial Point (m): 60196.82 60196.82 47062.33 46810.4 46810.4
##                                    X78      X79      X80      X81      X82
## Dist. from Initial Point (m): 45785.85 47834.95 47062.33 47062.33 47062.33
##                                   X83     X84     X85      X86      X87
## Dist. from Initial Point (m): 56837.6 56837.6 56837.6 55040.43 56064.98
##                                    X88      X89      X90      X91      X92
## Dist. from Initial Point (m): 55813.05 56064.98 56064.98 56064.98 56064.98
##                                    X93      X94      X95      X96      X97
## Dist. from Initial Point (m): 56064.98 56064.98 29292.15 29292.15 29292.15
##                                    X98      X99     X100     X101     X102
## Dist. from Initial Point (m): 55040.43 55040.43 55040.43 55040.43 55040.43
##                                   X103     X104     X105     X106     X107
## Dist. from Initial Point (m): 15689.43 15689.43 15016.56 15368.23 15368.23
##                                   X108     X109     X110     X111     X112
## Dist. from Initial Point (m): 45281.99 45030.06 45030.06 45802.67 45802.67
##                                   X113     X114     X115     X116     X117
## Dist. from Initial Point (m): 45030.06 46827.22 46827.22 46827.22 47079.15
##                                   X118     X119    X120     X121     X122
## Dist. from Initial Point (m): 45802.67 45802.67 46054.6 39302.52 39302.52
##                                   X123     X124     X125     X126     X127
## Dist. from Initial Point (m): 17215.86 16443.24 17467.79 16443.24 17467.79
##                                   X128     X129     X130     X131     X132
## Dist. from Initial Point (m): 17215.86 16695.17 16695.17 16695.17 16695.17
##                                   X133     X134     X135     X136     X137
## Dist. from Initial Point (m): 16695.17 16695.17 33608.63 33608.63 33608.63
##                                   X138     X139     X140     X141    X142
## Dist. from Initial Point (m): 33860.56 53512.02 53512.02 53763.95 44761.3
##                                  X143    X144     X145     X146    X147    X148
## Dist. from Initial Point (m): 44761.3 44761.3 44509.37 45785.85 44761.3 44761.3
##                                  X149     X150    X151     X152     X153
## Dist. from Initial Point (m): 44761.3 44509.37 10027.2 9775.264 9775.264
##                                   X154     X155    X156     X157     X158
## Dist. from Initial Point (m): 53008.15 53243.27 53495.2 47834.95 16695.17
##                                  X159     X160     X161     X162     X163
## Dist. from Initial Point (m): 16947.1 19264.95 18492.34 18492.34 18492.34
##                                   X164     X165     X166     X167     X168
## Dist. from Initial Point (m): 18492.34 19264.95 19264.95 18492.34 18492.34
##                                   X169     X170     X171     X172     X173
## Dist. from Initial Point (m): 17719.72 18492.34 19264.95 15385.05 15636.98
##                                   X174     X175     X176     X177     X178
## Dist. from Initial Point (m): 15133.12 15133.12 15133.12 15133.12 15133.12
##                                   X179     X180     X181     X182     X183
## Dist. from Initial Point (m): 15385.05 15385.05 55544.29 55544.29 55544.29
##                                   X184     X185     X186     X187     X188
## Dist. from Initial Point (m): 55544.29 55544.29 55292.36 55292.36 55292.36
##                                   X189     X190     X191     X192     X193
## Dist. from Initial Point (m): 56316.91 44240.62 18509.16 19533.71 19533.71
##                                   X194     X195     X196     X197    X198
## Dist. from Initial Point (m): 19533.71 18509.16 44811.77 44811.77 45063.7
##                                   X199     X200     X201     X202     X203
## Dist. from Initial Point (m): 44291.08 15658.97 16010.64 16683.51 16683.51
##                                   X204     X205     X206     X207     X208
## Dist. from Initial Point (m): 16683.51 16683.51 18059.73 18059.73 18059.73
##                                   X209     X210     X211     X212     X213
## Dist. from Initial Point (m): 18059.73 18059.73 13992.01 13992.01 39033.77
##                                   X214     X215     X216     X217     X218
## Dist. from Initial Point (m): 39033.77 39033.77 38781.84 38261.15 38261.15
##                                   X219     X220     X221     X222     X223
## Dist. from Initial Point (m): 38781.84 38781.84 11034.92 55040.43 55040.43
##                                   X224     X225     X226     X227     X228
## Dist. from Initial Point (m): 55040.43 55040.43 55040.43 12585.32 12585.32
##                                   X229     X230    X231    X232     X233
## Dist. from Initial Point (m): 12936.99 12585.32 34112.5 34112.5 33860.56
##                                   X234     X235     X236     X237     X238
## Dist. from Initial Point (m): 33087.95 25949.75 36195.24 36195.24 36195.24
##                                   X239    X240    X241    X242    X243     X244
## Dist. from Initial Point (m): 36195.24 35943.3 35943.3 35943.3 35943.3 46541.65
##                                   X245     X246     X247     X248     X249
## Dist. from Initial Point (m): 46541.65 46541.65 47314.26 50220.08 17938.01
##                                   X250     X251     X252     X253     X254
## Dist. from Initial Point (m): 17938.01 17938.01 17938.01 17938.01 17938.01
##                                   X255     X256     X257     X258     X259
## Dist. from Initial Point (m): 17938.01 17938.01 17988.47 17988.47 45584.38
##                                   X260     X261     X262     X263     X264
## Dist. from Initial Point (m): 45584.38 45836.32 45836.32 45315.63 45315.63
##                                   X265     X266     X267     X268     X269
## Dist. from Initial Point (m): 45315.63 46088.25 47381.55 48154.17 48154.17
##                                   X270     X271     X272    X273    X274
## Dist. from Initial Point (m): 48154.17 54519.75 39789.57 40041.5 40041.5
##                                  X275     X276     X277     X278     X279
## Dist. from Initial Point (m): 40041.5 40562.18 40562.18 40562.18 39789.57
##                                   X280     X281     X282     X283     X284
## Dist. from Initial Point (m): 39789.57 39789.57 39789.57 39789.57 39789.57
##                                   X285     X286    X287     X288     X289
## Dist. from Initial Point (m): 39789.57 38765.02 40041.5 51714.85 51714.85
##                                   X290     X291     X292     X293     X294
## Dist. from Initial Point (m): 51714.85 51714.85 59944.88 60196.82 59944.88
##                                   X295     X296     X297     X298     X299
## Dist. from Initial Point (m): 60196.82 10766.17 10766.17 10766.17 52773.04
##                                   X300     X301     X302     X303     X304
## Dist. from Initial Point (m): 53763.95 53763.95 54015.88 54015.88 54267.81
##                                   X305     X306     X307     X308     X309
## Dist. from Initial Point (m): 54015.88 54267.81 54015.88 54015.88 54015.88
##                                   X310    X311     X312     X313     X314
## Dist. from Initial Point (m): 53243.27 15653.8 13319.13 13319.13 13319.13
##                                  X315     X316     X317     X318     X319
## Dist. from Initial Point (m): 15653.8 14091.75 13839.82 14612.43 13319.13
##                                   X320     X321     X322     X323     X324
## Dist. from Initial Point (m): 45802.67 45802.67 45802.67 28771.46 28771.46
##                                   X325     X326    X327    X328     X329
## Dist. from Initial Point (m): 28519.53 53512.02 52739.4 52739.4 53512.02
##                                  X330    X331     X332     X333     X334
## Dist. from Initial Point (m): 52739.4 52739.4 53512.02 53512.02 53260.09
##                                   X335     X336     X337     X338     X339
## Dist. from Initial Point (m): 53512.02 16392.78 16392.78 16392.78 16392.78
##                                   X340     X341     X342     X343     X344
## Dist. from Initial Point (m): 16392.78 16392.78 16392.78 16392.78 17065.65
##                                   X345     X346     X347     X348     X349
## Dist. from Initial Point (m): 16392.78 56064.98 56316.91 57089.53 57089.53
##                                   X350     X351    X352    X353    X354    X355
## Dist. from Initial Point (m): 57089.53 57089.53 22326.8 22326.8 22326.8 22326.8
##                                  X356    X357    X358    X359    X360     X361
## Dist. from Initial Point (m): 22326.8 22326.8 22326.8 52739.4 52739.4 52991.33
##                                   X362
## Dist. from Initial Point (m): 52470.65

‘igraph’ package using the function ‘distances’ Create a hexagonal grid with 250m between hexagons

#create a hexagonal grid of Ada and Elmore Counties (takes about 20 minutes to complete)
#as you increase the resolution (decrease 'cellsize') the less duplicate verticies will be created
hex_studyCounties <- st_make_grid(studyCounties, cellsize = 250, square = FALSE)
plot(hex_studyCounties)

Create an ‘igraph’ network of neighbor cells

##create an 'igraph' network of neighboring cells which has a possiblity of being traveled through
#Steps to create an 'igraph' network
#1. create an edgelist - list to ‘from’ and ‘to’ cell id’s that represent neighbors 
edgeList <- as.matrix(as.data.frame(st_intersects(hex_studyCounties)))
#2. create graph of edgelist
edgeList_graph <- graph_from_edgelist(edgeList) #graph only shows the distance between cells (250 meters as defined above 'cellsize = 250')
#3. set the weight of all edges
E(edgeList_graph)$weight <- 250 #this could be adjusted by different factors such as elevation, habitat resistance, etc.

##run Least Cost Distance analysis on 'igraph' network
#orient hexagonal grid to the 'ABS' locals
orient <- as.numeric(st_intersects(coords, hex_studyCounties))
#remove any duplicate vertices (the higher the resolution the lower the number of duplicate vertices)
orient <- unique(x = orient)

Least Cost Distance calculation with ‘gdistance’

#run igraph LCD
studyCounties_leastCost2 <- distances(edgeList_graph, orient[1], orient[-1])

#put LCD calculation into data.frame
studyCounties_leastCost2_df <- data.frame(studyCounties_leastCost2, row.names = "Dist. from Initial Point (m):")
head(studyCounties_leastCost2_df)
##                                X1  X2  X3  X4  X5   X6   X7   X8   X9  X10
## Dist. from Initial Point (m): 250 250 500 500 750 1250 5250 5250 5500 5750
##                                 X11   X12   X13   X14   X15   X16   X17   X18
## Dist. from Initial Point (m): 60250 60000 59250 41250 41500 41250 58000 57750
##                                 X19   X20   X21   X22   X23   X24   X25   X26
## Dist. from Initial Point (m): 57500 57250 57000 56750 20250 20250 57500 57750
##                                 X27   X28   X29   X30   X31   X32   X33   X34
## Dist. from Initial Point (m): 57750 56750 58000 58000 57000 56500 57500 53500
##                                 X35   X36   X37   X38   X39   X40   X41   X42
## Dist. from Initial Point (m): 59000 56250 55750 54500 54000 53750 53500 53000
##                                 X43   X44   X45   X46   X47   X48   X49   X50
## Dist. from Initial Point (m): 52000 51750 51500 51500 28750 28500 28250 28000
##                                 X51   X52   X53   X54   X55   X56   X57   X58
## Dist. from Initial Point (m): 20000 20250 15500 66500 66250 66000 65750 50500
##                                 X59   X60   X61   X62   X63   X64   X65   X66
## Dist. from Initial Point (m): 50250 50250 50000 50500 61000 61250 61250 59250
##                                 X67   X68   X69   X70   X71   X72   X73   X74
## Dist. from Initial Point (m): 60750 60250 60500 30000 30000 30250 59750 59750
##                                 X75   X76   X77   X78   X79   X80   X81   X82
## Dist. from Initial Point (m): 17250 17250 17250 49750 49750 50000 50250 50500
##                                 X83   X84   X85   X86   X87   X88   X89   X90
## Dist. from Initial Point (m): 50000 51500 51500 51750 51750 50750 51000 40000
##                                 X91   X92   X93   X94   X95   X96   X97   X98
## Dist. from Initial Point (m): 40250 17750 17500 17750 17500 18000 17000 17250
##                                 X99  X100  X101  X102  X103  X104  X105  X106
## Dist. from Initial Point (m): 17250 32000 32000 32250 59000 49000 48750 48500
##                                X107  X108  X109  X110  X111  X112  X113  X114
## Dist. from Initial Point (m): 48750 49250 48500 49250 10750 10750 10500 58750
##                                X115  X116  X117  X118  X119  X120  X121  X122
## Dist. from Initial Point (m): 58000 58000 17000 17500 19500 19000 19250 19250
##                                X123  X124  X125  X126  X127  X128  X129  X130
## Dist. from Initial Point (m): 19500 18500 18750 18000 19250 15750 16000 15250
##                                X131  X132  X133  X134  X135  X136  X137  X138
## Dist. from Initial Point (m): 15000 14750 15750 60000 60000 61000 48000 19750
##                                X139  X140  X141  X142  X143  X144  X145  X146
## Dist. from Initial Point (m): 20000 51500 51500 51250 51000 17500 17750 18000
##                                X147  X148  X149  X150  X151  X152  X153  X154
## Dist. from Initial Point (m): 18250 19750 19750 19750 15500 15750 38750 38750
##                                X155  X156  X157  X158  X159  X160  X161  X162
## Dist. from Initial Point (m): 39000 38500 10750 60000 13250 13750 13750 13500
##                                X163  X164  X165  X166  X167  X168  X169  X170
## Dist. from Initial Point (m): 32500 32250 32000 31500 26500 35750 35500 35750
##                                X171  X172  X173  X174  X175  X176  X177  X178
## Dist. from Initial Point (m): 50000 49750 50500 57500 19000 19250 19250 19250
##                                X179  X180  X181  X182  X183  X184  X185  X186
## Dist. from Initial Point (m): 18750 19000 52250 52500 52750 52500 51500 52000
##                                X187  X188  X189  X190  X191  X192  X193  X194
## Dist. from Initial Point (m): 52000 54750 55000 55250 58750 39250 39250 39500
##                                X195  X196  X197  X198  X199  X200  X201  X202
## Dist. from Initial Point (m): 39250 39000 39000 38750 38750 38500 57000 56750
##                                X203  X204  X205  X206  X207  X208  X209  X210
## Dist. from Initial Point (m): 56500 56750 66000 65750 65500 11000 11000 59750
##                                X211  X212  X213  X214  X215  X216  X217  X218
## Dist. from Initial Point (m): 58250 59000 58750 59000 58750 58500 58250 15000
##                                X219  X220  X221  X222  X223  X224  X225  X226
## Dist. from Initial Point (m): 15250 15250 15250 15000 15000 15000 15250 29250
##                                X227  X228  X229  X230  X231  X232  X233  X234
## Dist. from Initial Point (m): 29500 58250 58250 58250 58000 58250 58500 58750
##                                X235  X236  X237  X238  X239  X240  X241  X242
## Dist. from Initial Point (m): 59000 58750 18500 18250 18500 18250 18500 18750
##                                X243  X244  X245  X246  X247  X248  X249  X250
## Dist. from Initial Point (m): 60250 60250 61750 23250 23500 23250 23500 57750
##                                X251
## Dist. from Initial Point (m): 56500

Benchmarking and Profiling - ‘igraph’ distances function & ‘gdistance’ costDistance fucntion

Benchmarking with the package ‘microbenchmark’

costDistance.time_mb <- microbenchmark(costDistance(studyCounties_ras_tr, fromCoords = as(as_Spatial(coords[1,]), "SpatialPoints"), toCoords = as(as_Spatial(coords[-1,]), "SpatialPoints")), times = 15)

distances.time_mb <- microbenchmark(distances(edgeList_graph, orient[1], orient[-1]), times = 15)

costDistance.time_mb
## Unit: milliseconds
##                                                                                                                                                               expr
##  costDistance(studyCounties_ras_tr, fromCoords = as(as_Spatial(coords[1,      ]), "SpatialPoints"), toCoords = as(as_Spatial(coords[-1,      ]), "SpatialPoints"))
##       min       lq     mean   median       uq      max neval
##  65.94644 66.31221 70.31113 66.84666 67.70158 99.20054    15
distances.time_mb
## Unit: milliseconds
##                                              expr      min       lq     mean
##  distances(edgeList_graph, orient[1], orient[-1]) 162.5931 179.4681 180.9126
##    median       uq     max neval
##  181.1937 183.1535 206.485    15
  • on average the function costDistance from the ‘gdistance’ package provided an output approximatly 100 milliseconds faster than the function distances from the ‘igraph’ package

Profiling with the package ‘profvis’

#perform profiling on 'gdistance' package
costDistance.prof <- profvis(expr = {
  #read in Idaho counties shape file
counties <- st_read("~/library/Mobile Documents/com~apple~CloudDocs/BSU/Graduate - MS Raptor biology/r_Projects/HES598/datatemp/original/idaho_cty.shp") %>% st_transform(crs = 32619)

#isolate Ada & Elmore countie pologons from Idaho county shape file
studyCounties <- subset(counties, NAME == "Ada" | NAME == "Elmore")

#bring in spatial coordiantes of 'ABS' locals
coords <- artBurrow_spdf

#set 'ABS' spatial coordinates from lat long to UTM 19N
coords <- st_as_sf(coords, coords = c("Longitude", "Latitude"), crs = 4326, stringAsFactors = FALSE) %>% st_transform(crs = 32619)

#create raster of 'studyCounties' spatial object
studyCounties_ras <- raster(x = extent(studyCounties), nrow = 250, ncol = 250)

#assign all cells to a value of 1 within the raster
studyCounties_ras <- rasterize(x = studyCounties, y = studyCounties_ras, field = 1)

#create an object of the Transition class
studyCounties_ras_tr <- transition(studyCounties_ras, transitionFunction = mean, 8)

#Geo correct transition object
studyCounties_ras_tr <- geoCorrection(studyCounties_ras_tr, type = "c")

#Conduct least cost distance calculation on all 'ABS' locals
studyCounties_leastCost <- costDistance(studyCounties_ras_tr, fromCoords = as(as_Spatial(coords[1,]), "SpatialPoints"), toCoords = as(as_Spatial(coords[-1,]), "SpatialPoints"))

}, interval = 0.01, prof_output = 'fstr-prof')
## Reading layer `idaho_cty' from data source `/Users/brentclark/Library/Mobile Documents/com~apple~CloudDocs/BSU/Graduate - MS Raptor biology/R_Projects/HES598/datatemp/original/idaho_cty.shp' using driver `ESRI Shapefile'
## Simple feature collection with 44 features and 17 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -117.243 ymin: 41.98818 xmax: -111.0435 ymax: 49.00091
## epsg (SRID):    NA
## proj4string:    +proj=longlat +ellps=GRS80 +no_defs
#perform profiling on 'igraph' package
distances.prof <- profvis(expr = {
  #create a hexagonal grid of Ada and Elmore Counties (takes about 20 minutes to complete)
#as you increase the resolution (decrease 'cellsize') the less duplicate verticies will be created
hex_studyCounties <- st_make_grid(studyCounties, cellsize = 250, square = FALSE)
#plot(hex_studyCounties)

##create an 'igraph' network of neighboring cells which has a possiblity of being traveled through
#Steps to create an 'igraph' network
#1. create an edgelist - list to ‘from’ and ‘to’ cell id’s that represent neighbors 
edgeList <- as.matrix(as.data.frame(st_intersects(hex_studyCounties)))
#2. create graph of edgelist
edgeList_graph <- graph_from_edgelist(edgeList) #graph only shows the distance between cells (250 meters as defined above 'cellsize = 250')
#3. set the weight of all edges
E(edgeList_graph)$weight <- 250 #this could be adjusted by different factors such as elevation, habitat resistance, etc.

##run Least Cost Distance analysis on 'igraph' network
#orient hexagonal grid to the 'ABS' locals
orient <- as.numeric(st_intersects(coords, hex_studyCounties))
#remove any duplicate vertices (the higher the resolution the lower the number of duplicate vertices)
orient <- unique(x = orient)

studyCounties_leastCost2 <- distances(edgeList_graph, orient[1], orient[-1])

}, interval = 0.01, prof_output = 'fstr-prof')

costDistance.prof
distances.prof

Introduce the packages

The critical step in the analysis is providing the least cost distances between all artificial burrow systems in the study area. Using the packages ‘igraph’ and ‘gdistance’ I can evaluate the least cost distance between each burrow geographic location utilizing their least cost distance fucntions (‘igraph’ - distances & ‘gdistance’ - costDistance).

R packages to try

“gistance” - Used to calculate distances and routes based on geographical grids. “igraph” - Used for network analysis on large graphs with millions of verticies and edges.

Outcomes

Both packages provided a comparable LCD’s for the ABS location data. The only thing of note is that when running the ‘igraph’ package it won’t allow duplicate verticies. These verticies must be removed before analysis can begin. It seems as though this can be mitigated by increasing the resolutoion (decreasing the cellsize) of the grid. Since many of the ABS locations are so close togother the lower resolution causes duplicate vertices/edges in the created hexagonal grid increasing the resolution of the grid will dramatically increase the amount of time that this analysis will take. For a dataset as small of this (>350 ABS locals) the increased processing time doesn’t seem appropriate.

Evaluate your choices

After benchmarking the ‘costDistance’ function from the ‘gdistance’ package the mean run time was approximately 80 - 130ms whereas as benchmarking the ‘distances’ function of the ‘igraph’ package provided a mean run time of approximately 200 - 230ms.

After profiling the ‘gdistance’ package and the ‘igraph’ package there was clearly a winner in efficiency. The ‘igraph’ package code script took nearly 280 times longer to run than the ‘gdistance package. The ’igraph’ package also used up more available memory than the ‘gdistance’ package due to the creation of a hexagonal grid which used approximately 300 Mb of memory and took the longest to run to completion.

There are definitely upsides to using both packages. The ‘gdistance’ package is ideal if the dataset (geographic coordinates) is relatively small whereas if the data being dealt with is very large then using the ‘igraph’ package would be more appropriate even though more time and memory allocaiton is required at the head of the project in order to create the graph.

As far as ase of use is concerned, the syntax of the ‘gdistance’ package is more intuitive (for me) and easy to use as there are far less functions available to sift through and use in the package as compared to the ‘igraph’ package. I learned how The amount of functions available for use in the ‘igraph’ packages reflects what the package is designed for, creating spatial networks. Due to the amount of functions avialable for the ‘igraph’ package it was difficult to keep track of what did what. For the purposes of my evaluation of the two packages, the ‘gdistance’ package makes more sense as the need for networking and the use of large datasets isn’t needed.

Reflect

Through this project I was able to learn more about landscape connectivity than I previoiusly knew before which will become more important for my future research of western burrowing owls in the state of Idaho. Even though I didn’t go into the complexities of least cost distance that come with an actual landscape, this exercise showed me how and helped me better understand what connectivity is and how it can be used. Before this assignment my experience and skill in the use of the ‘r’ programming language was very minimal and by working through the different packages I was able to better understand the basic syntax associated with ‘r’ and how to reseearch the different ways in which to find solutions to the problems that I faced.